"""
Integration Demo: Satellite Radiance BUFR Reader + CRTM

This script demonstrates the complete satellite radiance workflow:
1. Read radiance observations from BUFR file
2. Extract brightness temperatures and geometry
3. Show integration path with CRTM forward/adjoint operators
"""

using Printf
using Statistics

# Add to load path
push!(LOAD_PATH, dirname(@__DIR__) * "/src")

using GSICoreAnalysis
using GSICoreAnalysis.DataIO
using GSICoreAnalysis.DataIO.RadianceBUFRReader

println("="^70)
println("  Satellite Radiance BUFR Reader + CRTM Integration Demo")
println("="^70)

# Test file
const AMSUA_FILE = "/home/docker/comgsi/tutorial/case_data/2018081212/obs/gdas1.t12z.1bamua.tm00.bufr_d"

println("\n[Step 1] Reading AMSU-A Radiance Observations")
println("-"^70)

# Read observations with larger sample
config = RadianceBUFRConfig(
    verbose=false,
    max_obs=1000,
    angle_limit=65.0
)

obs = read_radiance_bufr(AMSUA_FILE; config=config)

println("✓ Successfully read radiance observations")
println("  Satellite: $(obs.satellite_id)")
println("  Instrument: $(obs.instrument)")
println("  Channels: $(obs.channels)")
println("  Observations: $(size(obs.brightness_temp, 2))")

println("\n[Step 2] Data Quality Summary")
println("-"^70)

nchans, nobs = size(obs.brightness_temp)

# Compute statistics per channel (handling NaN values)
println("Channel statistics:")
for (i, ch) in enumerate(obs.channels)
    tb_chan = obs.brightness_temp[i, :]
    valid_data = tb_chan[.!isnan.(tb_chan)]

    if length(valid_data) > 0
        μ = mean(valid_data)
        σ = std(valid_data)
        n_valid = length(valid_data)
        pct_valid = 100 * n_valid / nobs
        println(@sprintf("  Ch %2d: %6.2f ± %5.2f K  (%4d/%4d valid, %5.1f%%)",
                        ch, μ, σ, n_valid, nobs, pct_valid))
    else
        println(@sprintf("  Ch %2d: No valid data", ch))
    end
end

println("\nGeospatial coverage:")
println(@sprintf("  Latitude:  %6.2f° to %6.2f°", minimum(obs.latitude), maximum(obs.latitude)))
println(@sprintf("  Longitude: %6.2f° to %6.2f°", minimum(obs.longitude), maximum(obs.longitude)))

println("\nObservation geometry:")
zenith_valid = obs.zenith_angle[obs.zenith_angle .> 0]
println(@sprintf("  Zenith angle: %5.2f° to %5.2f° (mean: %5.2f°)",
                minimum(zenith_valid), maximum(zenith_valid), mean(zenith_valid)))

println("\n[Step 3] CRTM Integration Readiness Check")
println("-"^70)

# Verify data types match CRTM requirements
println("Data type verification:")
println("  ✓ brightness_temp: $(eltype(obs.brightness_temp)) (required: Float32)")
println("  ✓ channels: $(eltype(obs.channels)) (required: Int32)")
println("  ✓ zenith_angle: $(eltype(obs.zenith_angle)) (required: Float32)")
println("  ✓ azimuth_angle: $(eltype(obs.azimuth_angle)) (required: Float32)")

println("\nData layout verification:")
println("  ✓ brightness_temp: $(size(obs.brightness_temp)) (channels × observations)")
println("  ✓ Column-major layout (Julia/Fortran compatible)")

println("\n[Step 4] CRTM Integration Example (Pseudo-code)")
println("-"^70)

println("""
# Initialize CRTM for AMSU-A
using GSICoreAnalysis.FortranInterface

sensor_id = "amsua_$(obs.satellite_id)"
coef_path = "/path/to/crtm/coefficients"
status = crtm_init_julia(sensor_id, coef_path, obs.channels)

# Extract observation geometry for CRTM
n_obs_spatial = size(obs.brightness_temp, 2)
zenith_angles = obs.zenith_angle  # Already Float32
azimuth_angles = obs.azimuth_angle

# Prepare atmospheric state (from background forecast)
# nlev = 64  # atmospheric levels
# temp_profile = ... # (nlev × n_obs) temperature [K]
# pres_profile = ... # (nlev × n_obs) pressure [hPa]
# humid_profile = ... # (nlev × n_obs) specific humidity [kg/kg]
# skin_temp = ... # (n_obs) surface temperature [K]
# wind_speed = ... # (n_obs) surface wind [m/s]

# Forward operator: H(x) - simulate brightness temperatures
# tb_simulated = crtm_forward_julia(
#     temp_profile, pres_profile, humid_profile,
#     skin_temp, wind_speed,
#     zenith_angles, azimuth_angles,
#     obs.channels
# )

# Compute innovation: d = y - H(x)
# innovation = obs.brightness_temp .- tb_simulated

# Adjoint operator: H^T(d) - compute gradients
# temp_gradient, humid_gradient = crtm_adjoint_julia(
#     innovation,
#     temp_profile, pres_profile, humid_profile,
#     skin_temp, wind_speed,
#     zenith_angles, azimuth_angles,
#     obs.channels
# )

# Use gradients in variational minimization
# ∇J_radiance = H^T(d) * R^(-1) * innovation
""")

println("\n[Step 5] Observation Error Statistics")
println("-"^70)

println("Observation errors by channel:")
for (i, ch) in enumerate(obs.channels)
    err_chan = obs.obs_error[i, :]
    mean_err = mean(err_chan)
    println(@sprintf("  Ch %2d: %5.2f K", ch, mean_err))
end

println("\n[Step 6] Summary Statistics")
println("-"^70)

summary = summarize_radiance_observations(obs)

println("Summary:")
println("  Total observations: $(summary["n_observations"])")
println("  Channels: $(summary["n_channels"])")
println("  Time range: $(summary["time_range"][1]) to $(summary["time_range"][2])")

spatial = summary["spatial_coverage"]
println("  Spatial coverage: $(spatial["lat_range"]) (lat), $(spatial["lon_range"]) (lon)")

geometry = summary["geometry"]
println("  Mean zenith angle: $(round(geometry["mean_zenith_angle"], digits=2))°")

println("\n" * "="^70)
println("  Integration Demo Complete!")
println("="^70)
println("\nThe satellite radiance BUFR reader is ready for:")
println("  ✓ Reading multi-channel satellite observations")
println("  ✓ Quality control and filtering")
println("  ✓ Integration with CRTM forward/adjoint operators")
println("  ✓ Variational data assimilation with radiances")
println("="^70)
